/****************************************************************************
 * Elisa 1.0                                                                *
 *--------------------------------------------------------------------------*
 * Author: Laurent Granvilliers, Vincent Sorin                              *
 * Copyright (c) 2004 Laboratoire d'Informatique de Nantes Atlantique       *
 *--------------------------------------------------------------------------*
 * Elisa is distributed WITHOUT ANY WARRANTY. Read the associated           *
 * COPYRIGHT file for more details.                                         *
 *--------------------------------------------------------------------------*
 * elisa_newton_reduction.h                                                 *
 ****************************************************************************/

/**
 * @file elisa_newton_reduction.h
 * @brief Implementation of the multidimensional interval Newton method.
 */


#ifndef _ELISA_NEWTON_REDUCTION_H
#define _ELISA_NEWTON_REDUCTION_H

#include "elisa_alloc.h"
#include <vector>
#include "elisa_defs.h"
#include "elisa_reduction.h"

#include "elisa_observer.h"

#include "elisa_real_constraint.h"
#include "elisa_matrix.h"
#include "elisa_interval.h"
#include "elisa_interval_test.h"
#include "elisa_interval_form.h"

#ifdef ELISA_PROFILE
#include "elisa_time.h"
#endif // ELISA_PROFILE


NAMESPACE_ELISA
/**
 * @ingroup Reduction
 * @class NewtonReduction
 * @brief Multidimensional interval Newton method
 */
class NewtonReduction : public Reduction {
public:
	/// Constructor
	NewtonReduction (VarImplSet & vs,
	                 SetConstraints & s,
	                 IntervalFormFactory *createForm);

	/// Destructor
	~NewtonReduction();

	/// Base class's virtual methods
	//@{
	double priority() const
	{
		return 10.0;
	}

	bool reductionAlgorithm();
	//@}

	/// Set the improvement factor
	void setImprovement(double p)
	{
		m_improve = p;
	}

	/// Set the improvement factor for all further created reductions
	static void SetImprovement(double p)
	{
		StaticImprove = p;
	}

#ifdef ELISA_PROFILE
	/**
 	 * Accumulates the time spent in the reduce operation
 	 * Only available in profiling mode.
 	 */
	static Clock ClockNewton;

	/**
 	 * Counts the number of calls of the reduce operation
 	 * Only available in profiling mode.
 	 */
	static long NumCall;
#endif



	string getReductionName() const
	{
		return "Newton";
	}


	void print(std::ostream & os, std::string indent, VecConstraints *v) const
	{
		os << "NewtonReduction [";
		printConstraints(os, "", v);
		os << "]";
	}

private:
	std::vector<IntervalForm *> _exprs; // interval form of F, given F=0
	std::vector<RealVarImpl *>   _vars; // The variables
	int _dim;                         // dimension of system

	double m_improve;               // improvement factor
	static double StaticImprove;  // default improvement factor

	// Particular matrices created only once
	DenseMatrix<Interval> _MIzero; // zero matrix of intervals
	DenseMatrix<Interval> _VIzero; // zero vector of intervals
	DenseMatrix<double>   _MRid;  // identity matrix of reals

	/// @name Auxiliary functions
	//@{

	/**
 	 * mid := interval midpoint of variable domains
 	 * b := -F(mid)
 	 */
	void evalMidpoint(DenseMatrix<Interval> & mid,
	                  DenseMatrix<Interval> & b) const;

	/// y := (variable domains) - mid
	void evalVariable(DenseMatrix<Interval> & y,
	                  DenseMatrix<Interval> & mid) const;

	/// J := Jacobian matrix
	void evalJacobian(DenseMatrix<Interval> & J) const;

	/// mid := midpoint of J
	void evalMidJacobian(DenseMatrix<double> & mid,
	                     DenseMatrix<Interval> & J) const;

	/// (variable domains) := y + mid
	bool reduceDomains(DenseMatrix<Interval> & y,
	                   DenseMatrix<Interval> & mid) const;

	/// One Gauss-Seidel iteration over the system Ax = b
	static bool GaussSeidel(DenseMatrix<Interval> & A,
	                        DenseMatrix<Interval> & x,
	                        DenseMatrix<Interval> & b);

	/// p := m-1
	static bool Inversion(DenseMatrix<double> & p,
	                      DenseMatrix<double> & m,
	                      DenseMatrix<double> & id);

	/// m := p * a
	static void Multiply(DenseMatrix<Interval> & m,
	                     DenseMatrix<double> & p,
	                     DenseMatrix<Interval> & a,
	                     DenseMatrix<Interval> & zero);

	/// Returns true if A is diagonally dominant
	static bool IsDiagonallyDominant(DenseMatrix<Interval> & A);
	//@}

	/// Copy protection
	NewtonReduction & operator=(const NewtonReduction & o)
	{
	}
};
END_NAMESPACE_ELISA

#endif // _ELISA_NEWTON_REDUCTION_H
